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A method is given for developing approximations to neutron and gamma-ray transport distributions 
in the form of a superposition of £/* polynomial series multiplied by exponentials and characterized 
by different scale factors. Convergence conditions and error bound expressions are given. Biorthogo- 
nality properties are worked out. Examples are given in which these approximations are compared 
with polynomial approximations. 

Key words: Biorthogonal function; gamma-ray transport; moment methods; neutron transport; poly- 
nomial representations; radiation penetration. 

1. Introduction 

In the moment method of calculating space distributions in radiation transport theory, spatial 
moments denned by 

1 f x 

M n = —\ dzz n F(z), (1) 

n ! Jo 

or 

M n = —^ ( X dzz»F(z), (V) 

are evaluated by numerical integration of a recursive system of Volterra integral equations of the 
second kind [1, 2]. 1 For any of several reasons, one frequently encounters the problem of construct- 
ing distributions by use of alternate, rather than consecutive moments [2,5]. Because spatial trends ■ 
are basically exponential, this leads to the possibility of representations in a family of biorthonormal 
systems of polynomials U k (z), with adjoint polynomials £/J[(z), satisfying the condition 



r 

n N ' m v 
Jo 



dz*er*U* a (z)U'Jx) = 8„„ (2) 

Jo 

with 

z^(z) = i (-)«(-) j^. (2<) 

Many other properties of these functions have been given, in references [1~4], and particularly 
in reference [4]. Perhaps we should say here that we draw freely on the results and terminology 
of reference [4] and reference [6] in this paper. 
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The application of these polynomials to gamma-ray and neutron distributions leads to questions 
(a) of convergence in the limit, and (b) of convergence of a truncated series to within known error 
bounds. In the latter case, the utility of the method depends on the rapidity of convergence, because 
it is not practical to evaluate more than 6 to 10 of the series coefficients. 

The problem of convergence in the limit was treated in reference [4], and a more recent 
paper [6] developed error bounds for truncated series. In addition, for gamma-ray distributions, 
the convergence does tend to be rapid. Thus the problem of constructing distributions from mo- 
ments would appear to be completely solved. The convergence arguments are summarized in 
appendix B, with additions. 

Unfortunately, however, many important problems do not have rapidly converging series. For 
example, if the radiation source has components of differing energy, their attenuation is a super- 
position of exponentials, and the more rapidly attenuated components give slowly converging 
contributions to the coefficients of the polynomial representation. 

Another case is that in which the radiation initially travels at an oblique angle to the direction 
of penetration. A still more important case is that of the neutrons, where very high numbers of 
collisions with nuclei give spatial distributions which combine exponential with Gaussian features 
[7]. 

Existence of these problems led very early to use of nonlinear approximations which were 
designated "function-fitting" methods, in which the individual terms in an approximating series 
have the same functional form, with differing weight and differing scale factor [2]. While very 
useful, and in good agreement with the more basic polynomial representation, these methods have 
the limitation that it is difficult to assess their accuracy. In fact, it is easy to show that the more 
obvious error bounds which can be determined for such representations must be broader than error 
bounds for the polynomial representation. This is highly unsatisfactory, because one turns to these 
methods mainly when the polynomial representations are poor. 

In this paper we develop a type of approximation which we designate a "plural-series" repre- 
sentation or approximation. It is a generalization of the polynomial approach which incorporates 
the flexibility of the function-fitting methods. But unlike the latter, one can prove convergence in 
the limit, and can evaluate bounds to truncation errors. Thus it combines the advantages of both 
methods. Because the polynomial method is a special case, we will refer to it herein as a "single- 
series" approximation. 

We will show that it is possible to use plural-series calculations to evaluate error bounds 
for at least some function-fitting representations. This makes it possible to evaluate accuracy 
for previously published calculations. 

We have begun to apply the plural-series method to types of distribution for which poly- 
nomial approximations converge poorly, and we include an example of this in which the distri- 
bution has no physical significance. 

2. Plural-Series Representations 

The representation which we wish to use in obtaining an approximation to F(z) is 

F K (z) = 2(Ajlx*")F{(z), (3) 

where 



F{{z) = ^ a{W n (zlxj)e-*l*j. (3 ') 



The Xj are arbitrarily chosen scale factors and the Aj are arbitrary multipliers. The approximate 
function Fn(z) is required to have TV moments which have values exactly agreeing with known 
values A/*, Affc+2, M^+4, . . .for the corresponding sequence of moments of F(z). 
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For J > 1, there are J - N coefficients a,J in (3), a much largej number than the number of 
moments, N. Hence other requirements must be invoked fo determine the coefficients. For this 
purpose we use a simple variational rule, that the sum S\ defined as follows is to be minimized: 



S N =% {AJxf+i) 2 (c n a{)\ 



71=0 



(4) 



where the c n are arbitrary weighting coefficients. This leads to an elementary Lagrange multiplier 
calculation for the coefficients a J w , with side conditions as follows expressing the reauirements 
of the moments of the representation: 



M Mn =2 (Ajlx^)M{ +2n , 



j=i 



vhere 



1 f°° 

MJ = dzz k + 2n FJ (z) 

w fc+2 „ (t + 2n)\Jo n() ' 



(*• 



(5) 



(5') 



Because of the simphcity of (2'), the expression for a J n in terms of the moments of the "partial' 
functions F J N is very simple. Since 



al = £° d(zlxj) (zlx } )*U*(zlxj)Fi N (z) , 



(6) 



insertion of the expression (2') for [/$ into (6) gives 



«i"2H'(J)«U'. 



A+2J+1 



(7) 



Before proceeding with the further analysis, let us mention restrictions to be observed on the 
various types of arbitrary coefficients. The Xj must be positive or else the integrals do not converge. 
In addition, we will limit our considerations to positive values oiAj so thatS.v will always be positive. 
Then Sn turns out to be a squared norm and plays a central role in the discussion of error bounds. 
In addition, we assign positive values to the weighting coefficients c n , which play a role here similar 
to that of the c n coefficients in reference [6]. 

What the minimizing condition on Sn does is to select a particular biorthogonal system for the 
approximation. Our task is largely that of familiarizing ourselves with properties of this system, 
which turn out to be attractive as well as useful. 

3. Evaluation of Coefficients 

Equation (7) expresses a linear system whose formal solution makes it easier to express the 
main calculation to be performed. In matrix form, (7) is seen to be simply 




/ ] 




1 -1 
1 -2 1 



1 



:IV 








x~r 

j 



I w " \ 

ML* 

Mi 
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Since the binomial matrix is itf own 









1 
x) 
x 4 



nverse, we can immediately write 


-1 

-2 1 



it\ 



a i 



(8) 



Designating the matrix elements of the product of the first two matrices on the right by y? , we 
have the desired result, 

N-l 

Ml =x k+1 y yi. aJ (9) 



Let us note in passing that the y j . are the coefficients which, according to (2'), change a 
series of E/J polynomials into a power series, 



,2i N-l 



I WW*). 



(10) 



(k +20! 

To obtain the form of the side conditions which we wish to use, we insert (9) into (5), with the 
suit 



j=l w=0 



(11) 



Then, designating the N Lagrange multipliers as k u X 2 , . . . h N , we write out the quantity to be 
minimized, 

S N = ± (Ajlxf+i) ff (<**)«}+ 2 2 X, (M, +2i - 2 4 J 7i«il • (12) 

j=l ^n=0 J i=0 ^ j=l n=0 ^ 

Vanishing of linear variations of Sn with respect to the Xj gives back equations (11). The conditions 
{dSN/da j n = 0} yield the following additional J • jV equations: 



N-l 

I 

i=0 



«i^rs^#n 



(13) 



Formal solution of the linear system of J • (A^-h 1) equations requires only a straightforward 
matrix inversion. We insert (13) into (11) to obtain 



N-l 

2 

i'=0 



Mk+2i — 2 Udi'^i'i 



where 



J N-l 

iHv = y Ajx k + l y yj y{ /c 2 . 

j=l w=0 



(14) 



(15) 
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The matrix whose elements are uu> is symmetric and readily evaluated. Writing v nn ' for the elements 
of the inverse matrix, the equations (14) become 

yv-i 

K= ^ Vnn'M k+2 n', (16) 

n' = 

and values of the a{ follow from (13). 

4. Norm Convergence for N— » °° 

Let us assume that 7=1, that x\ and c n values have been selected, and that for this case 
lim Sn is finite. Then it follows that for J > 1, for the same set of c n values, and for a set of xj values 

which include the same Xi, the sequence of Sn values must again have a finite limit. In addition, 
the sequence of Sn values is monotone increasing. 

To demonstrate that these assertions hold true we first note that Fn+i is an approximation 
which fits the same moments as does Fn (plus one more). If, in the expression for F.\+i, one arbi- 
trarily changes the a\, values to zero, for ally, the first /V moments are unaffected due to biortho- 
gonality. Call such a modified approximation F N +i, and calculate the corresponding quantity 
S'n+i' Then 

Sn + i^S n+ i^S n , J 2*1. (17) 

The first of the inequalities holds because the aL would not usually be zero, while the second 
holds because S.v was calculated to have the smallest value possible for a set of approximations 
which includes Sjv+i. Thus the S\ sequence is monotone increasing, as well as positive. 

That the Sn sequence for J > 1 is bounded follows if the J— I sequence converges, and if 
any one of the Xj has the same value in the plural series case as the xi of the single-series case, as 
we have assumed. After all, for each value of /V, the single-series approximation is only one in 
which a-* ' = for ally except that corresponding to the single-series X\. It is therefore an element 
in the set of approximations among which the plural-series Sn is smallest. Hence the value of Sn, 
J > 1, is always bounded above by the single-series limit and therefore the limit, TV— » °°, exists 
and is less than or equal to the single-series limit. 

This argument makes no requirement on the values of Xj, J > 1, other than convergence of a 
single-series sequence in which one of the Xj values is used as scale factor. Therefore, plural- 
series approximations can include components with exponential trends inappropriate for use in 
single-series representations. 

5. Transformation Properties 

The transformation properties of plural-series approximations follow from the fact that adjoint 
to Fa (z) is a function which we designate F A , 

We establish this by an argument in terms of "transformation kernels", and for this purpose we 
first observe that 

{ci»£/*(z)e-} = ^ dyy k e<*+»T*{z,y){cnU*i{y)}, (19) 
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where 

T k (z,y) = 2 ic-lUi(z)}{c-lUXy)}. (19') 

tt=0 

The series (19') is assumed either to converge for all y, z or else to diverge at z=0, y=0 in such 
a way that the integration (19) still leads to convergent results. For c n =l, T k is a combination 
of modified Bessel functions, with a divergence at y = z = [4]. Because of the factor y k , this diver- 
gence is harmless in the integration. 

If we insert (10) into (18), and then mdke use of (13), we obtain the expression 

F K (z) = */<*"> ^ < c » V#*l*i) ■ (20) 

From (19) it then follows that 

FUz) = ( X dyy k e-^y^jT k (zlxj,ylxj)F N (z), (21) 

Jo 

so that 

F N (z) = j* dyy* { £ (A J l**+ 1 )e-<*+*l*jl*(zlx j9 ylxj)] F N (y). (22) 

Thus, the transformation appropriate to the plural-series representation has a kernel which is 
a simple linear combination of the kernels for the component series taken separately. 

To conclude this argument, we note that multiplication of (13) by (c«) 2 a J , followed by sum- 
mation over n, yields 

2 (cna{ ) 2 = 2 XiMU = f X dz2*F N {z)F{ (z) , (23) 

rt=0 i=0 ^° 

S N = r dzz k F N (z)F N (z). (24) 

Jo 



that 



6. Bio'thogonality Properties of Plural-Series Systems 

The discussion of error bounds is perhaps simplest in terms of the biorthogonal functions 
appropriate to plural-series cases. Let us therefore derive expressions for the two sets of biorthog- 
onal functions, which we designate P.x(z) and Py(z). 

We consider first the sequence of moments Mk= 1, and Mk+2n =z 0, for n > 0. The "function" 
which has these moments is irrelevant, though it must basically be a Dirac delta function located 
at z=0. In any case, when this "function" is represented in any orthogonal or biorthogonal sys- 
tem, we expect all the terms to appear. The sequence of approximate representations correspond- 
ing to this moment system will be designated Gn(z); and a corresponding sequence of adjoint 
functions, Gn(z) , can be constructed by use of (18). 

Two sequences of difference functions give the biorthogonal functions which we seek. For 
Nz* 1, 

Pn(z) = {G x+1 - G N }ld N , (25) 

and 

Pn(z) = {G N+1 - G N }/d N , (26) 
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where ds is a normalizing constant still to be specified. 

Since Cjv+i and Gn have identical moments Mk+2n, for n = 0, 2, . . . 2(N—\), it is clear that 
Pn is orthogonal to P n for n < N, the latter being a linear combination of even powers: 



Jo 



dzz k P n (z)Ps(z) =0,n< N. (27) 

To prove that P\(z) is also orthogonal to P n , for n < /V, we write 

T dzz k P s Pn = \"dzz k p dyy k {^ (AjIx^THz/xjMe^M'jl Px(z)P n (y) 



j=t 



Jo 



dyy k P n {y)Px(y) = 0. (27') 



Therefore (25) and (26) describe a biorthogonal system which obeys the transformation (22). To 
complete the system we add the terms 

Po(z) =Gi(z)ldo, (28) 

Po(z) =G 1 (z)ld . (28') 

The normalizing condition for this system is 



f 

Jo 



dzz k P s {z)Ps(z) = \, (29) 

which may be written by means of (24) and (4) in the form 

d l= 2 ^P S <#(«*>« - KM', (30) 

j=l n=0 

.7=1 

where («{).v is the coefficient a£, determined for the construction G#, and (a J v ).v is always zero. 
This completes the specification of the plural series biorthonormal system. In the special 
single-series case, there is cancellation between the quantities in braces in (30); and the usual 
normalization is recovered if we choose c n = 1, and Ai=x k+l . 

7. Convergence of Plural-Series Approximations 

Norm convergence of plural-series approximations was earlier shown to hold if there is norm 
convergence of any of the single-series approximations which can be formed using the plural- 
series scale factors Xj. But this does not guarantee pointwise convergence of these approxima- 
tions. To examine the latter, let us first note that (3) may be rewritten in the form of a Pn series, 



F. v =2/„P„(z), (31) 

/» = () 

with 
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Il^.v|| 2 = 2 f* = S »- (32) 

H = 

The form (31) follows from (22); one must recall that F\ is a polynomial and can therefore be 
written as a linear combination of the polynomials P n . Equation (32) follows from the normaliza- 
tion (29), 

||P.v|| 2 =l. (33) 

We proceed by designating the kernel function of the transformation (22) as $ rk (z, y), with 
the following extension of the notation to the A^'th finite sum: 

F%{z, y) = 2 (Ajlx^ 1 ) J c^U k n (zlxj)U k n (ylxj)e-^^ )lx K (34) 

j=l n=0 

Next, applying (22) to P n , but noting that U n '(zlxj) is orthogonal to the polynomial P n (z) for 
ally and all n' ^ N > n, we observe that 

P n (y)= r dzz k 3 rk (z,y)P n (z)= r dzz*lr k N {z,y)P n {z). (35) 

Jo Jo 

The second equality holds for n < N, and due to biorthonormality, 3^m(z, y) must have the form 

WHz,y)=^Pn(y)Pn(z) + '2 r n (y)P n (z), (36) 



where we do not have to concern ourselves with values for the coefficients r n {y) of the remainder 
term. From (36) and (32) it follows that for fixed y, 

2[ /> -<y)] , *ll^&ll 1< ^'*0 r .y). (37) 

since ST k {y, y) is identical with || ST k \. 

Using this result together with the Schwarz inequality, for z ¥" 0, we see that if the set of 
coefficients f n in (31) and (32) has the property Sn ^ S x < », then 



n=M ^ n=M J n=M rc=M 



(38) 



The term on the right goes to zero as M, TV-* o°; hence the Cauchy criterion is satisfied, and 
S/w^w converges. Thus, pointwise convergence of a plural-series approximation follows if there 
is norm convergence of any of the single-series representations based on the plural-series scale 
factors %j. 

There remains the question as to whether the limit of a plural-series sequence of approxima- 
tions agrees with that of the single-series approximations. To establish this we refer to the form 
(3), (3'), obtainable for any representation (31). We first note that \i xj is the largest of the plural- 
series scale factors, then the functions 

Fi{z) = 2 a{U>< n (zlxj)e-*l*Jj=l,2, . . ., 7, (39) 

n=0 

in any convergent representation of the type (3), (3'), can also be represented by a convergent 
series using C/J (z/xj) only. One has, therefore, a convergent single-series representation for 
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2 (AJx^FHz), (39') 

if one exists for F(z). The coefficients of the single-series representation of this function must 
be exactly those of the corresponding single-series representation for F(z). because both are 
given by the same linear combinations of the moments of F(z) if the same series representation 
(based on xj) is used for both. Hence (39') is equivalent to a single-series representation of F(z) 
in the limit. 2 

8. Bounds to Truncation Errors 

To calculate bounds to truncation errors we require norm convergence (37), which follows 
from norm convergence of a single-series representation based on one of the xj 9 as already stated. 
As a result of (35), the kernel of the transformation (22) must be representable in a form analogous 
to (19'), 

^~ k (z,y)= 2 (Ajlx^)e^^hT k (zlx h ylxj) = j? P n (z)P n (y). (40) 

j=l «=0 

This follows from biorthogonality, plus the fact that, considered as a function of z, say, for fixed 
y, this function can be represented as a series of U k n (zlxj) functions with finite norm, where xj 
is the largest scale factor. (See appendix B.) 

Following the notation of reference [6], we put y—z and make the designation 

[Dl(z)Y = ^ (Ajlx^)e-^j THz/xj, zlxj) = f [P„(z)] 2 . (41) 

j=l 11=0 

For approximations correct to N terms, one can accurately use the following truncated kernel 
function: 

Y P„{z)P„(y). 

SA (42) 

If we now apply the Schwartz inequality to the quantity (F — F„) as in (19') of reference [6], 
we obtain 



^ ». — ft J « — V ^ » — ft J 



(43) 



or 



\¥-F*\* VsZD$(z), (43') 



where we employ the notation of reference [6] except that exponentials now appear both in the 
function being approximated and in the definition of Z)^(z), 

[£>*(;)]> = [D*(z)]» - J [Pn(z)Y. (44) 

n=0 

We can follow reference [6] farther by defining an "error build-up" function Eb(z) , 



1 See Appendix B for a discussion of the automatic exclusion of functions having all moments zero. 
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Eb(z) = VS x D*(z)IF(z), 
so that the (fractional) bounds to the truncation error are given by 



\F-Fs\ 



= £6(2) 






(45) 



(46) 



Both factors on the right of (46) reduce in the single-series case to the results of reference [6] 
because the extra exponential factor cancels out. 

Examples: In transport problems, different energy components, characterized by exponentials 
which decrease at different rates, are commonly superposed. We choose a very simple case of 
this type as an example, namely 



g(z) = 0.75e-*/°- 2 + 0.25e~*. 



(47) 



Table 1 gives the function e z g(z), together with five approximations, of which (a) is a single-series 
case, while (b) is a 3-term plural-series case. The exponential coefficients Xj are listed below the 
table of values. 3 The plural-series approximation (b) has a component with one scale factor very 
near unity and another which is about 0.33, and therefore much closer to the value 0.2 in (47). 
All calculations were performed with the computer program sketched in appendix A. 

Below the scale factors in table 1 and norms Sn are listed for 5 through 8 terms in the approxi- 



3 See Appendix A for the rules used in selecting both xj and Aj values. 

Table 1. Values of g(z)e z and five approximations. 
Scale factors Xj are listed below the values of the function. At the bottom of the table are the last four values of Sy together with an extrapolated value for S x 







(a) 


(b) 


(c) 


(d) 


(e) 


z 


Exact 


k = 


k = 


k = 2 


£=6 


k = 6 






7=1 


7=3 


7 = 3 


7=3 


7=3 





1.0000 


0.680 


0.799 


1.000 


1.000 


1.000 


0.05 


0.8641 


.650 


.737 


0.864 


0.863 


0.864 


0.10 


.7527 


.622 


.682 


.753 


.748 


.753 


0.2 


.5870 


.569 


.588 


.587 


.574 


.589 


0.4 


.4014 


.478 


.452 


.402 


.377 


.406 


0.6 


.3180 


.406 


.366 


.318 


.292 


.323 


0.8 


.2806 


.349 


.312 


.280 


.252 


.284 


1.0 


.2637 


.305 


.278 


.263 


.246 


.265 


1.5 


.2519 


.241 


.242 


.252 


.243 


.249 


2.0 


.2503 


.220 


.237 


.251 


.248 


.251 


3. 


.2500 


.238 


.248 


.250 


.254 


.249 


5. 


.2500 


.271 


.256 


.251 


.248 


.250 


7. 


.2500 


.232 


.243 


.250 


.251 


.250 


10. 


.2500 


.255 


.254 


.251 


.250 


.250 


13. 


.2500 


.324 


.267 


.253 


.249 


.250 


16. 


.2500 


.0835 


.199 


.241 


.254 


.250 


20. 


.2500 


.0717 


.241 


.264 


.243 


.249 


24. 


.2500 


3.505 


1.08 


.297 


.251 


.252 


28. 


.2500 






.722 


.297 


.254 


34. 


.2500 






-1.06 


.077 


.226 




X\ 


1.000 


0.991 


0.801 


0.850 


0.771 




Xi 




.572 


.462 


.491 


.474 




X3 




.330 


.267 


.283 


.280 




5*, #=5 


0.234 


0.2761 


0.7506 


0.7606 


0.01013 




6 


.249 


.2789 


.7514 


.7611 


.01014 




7 


.262 


.2821 


.7521 


.7611 


.01023 




8 


.275 


.2856 


.7524 


.7613 


.01023 


•j 00, 


extrapolated: 


(0.446) 


(0.324) 


(0.759) 


(0.763) 


(0.0108) 
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mation. The final value in parentheses represents a very rough (linear) extrapolation to S^ using 
yy-1/2 as variable; while arbitrary, it tends to be quite conservative: 



c _ 

■JN+l 



V/V+l 



Sn-i 



v/v+i v/v-i. 



(48) 



The difference between the 8th norm and the extrapolation is an indicator of the degree of con- 
vergence, although the extrapolated value cannot be taken very seriously, particularly when, as 
in case (b), successive differences are still increasing. It turned out that with increasing iV, these 
extrapolations first increased, and then in high-quality approximations decreased to meet the 
(rising) norm values. 

Both the approximation values and the data on norms indicate that the plural-series approxi- 
mation is better, the single-series approximation being generally unacceptable. 

The remaining three columns represent an attempt to improve the approximations through 
use of information about the value and the first derivative of g(z) at z = 0. This information makes 
possible the approximation of a modified function, selected to be 

g'( z )= g ( z )-e-*i^ (49) 

which has a trend near z=0 which is proportional to z 2 , and which can thus be represented with 
plural-series approximations involving U% polynomials with k ^ 2. Three-term plural-series approxi- 
mations are used in all cases (c)-(e). The cases (d) and (e) ignore the first four even moments and 
use f/g polynomials [9]. In addition, for case (e) the Xj values were scaled down according to the rule 



=,a = 0.3. (50) 

j 

This form is suggested by a condition that a single adjoint function of the form cos(az) be trans- 
formed by the appropriate kernel functions T(z/xj, y/x'j), as in (21), to the exponentials e~ zlxj . 

Some of the sharp behavior for z near zero is subtracted out by the additional term in (49), so 
that one expects improvement in the convergence, in the values for large z, and especially for small 
z which is now given correctly in the limit. But because there is sharp behavior remaining, one does 
not expect the improvement to be as dramatic as shown. 

We should note that reducing the scale factor would also improve (a) and (6); the advantages in 
rising build-up factors which were discussed in references [6] and [8] apply here also. 

Dropping the lower moments tends to degrade the approximation (d) relative to the approxima- 
tion (c) for small z values, while introducing an improvement for large z values. Reducing the Xj 
values, on the other hand, improves the quality at all z. Recall that at least one of the Xj must always 
give a converging single-series approximation. That is the case for all approximations here shown. 

Table 2 compares errors with estimated error bounds. 4 In all cases shown, we have used Sg 
rather than the extrapolated approximation to S x as multiplier in estimating error bounds, ^s a 
result, in (a) at z = 0A the estimated error bound falls slightly below the error. This shows the 
problem of "estimating" an error bound in a poorly convergent case. Perhaps the extrapolated 
Soo would be a better multiplier to use generally, though there is not much difference when con- 
vergence is at all satisfactory. 

Elsewhere the estimated error bounds are reasonable, even excessive. Some use of the option 
c n # 1 was made, with the following prescription [6]: 

c„-(/i+l) m/4 , (51) 



4 See appendix A for details on the calculations. They made use of (44) to (46), with data on D? obtained using(19') with y = z. 
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Table 2. Comparisons between the fractional error, (F — F N )/F, designated "Err," and an estimated fractional error bound, 

designated "EB." 
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-.049 
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with m = for cases (a) and (6), and m= 2, 3, 6 for cases (c), (d), and (e), respectively. Large values 
of m can lead to poor convergence of the norms, but for the functions g(z) and g ' (z) approximated 
here, choice of m does not affect ultimate convergence. The evidence of table 1 is that the choices 
were reasonable. 

Unlike the single-series case, the value of m affects the approximation as well as the error bound 
functions. 

Error bounds generally tend to become large as one approaches large values of z at which 
convergence breaks down. But one also finds extremely large error bounds for small values of z 
in cases (d) and (e). These are due to the insensitivity of large moments to the details of the distri- 
bution at small z. That these bounds are realistic was demonstrated in some calculations not 
shown, in which for some m values, with calculations otherwise like (d) and (e), the errors in the 
small z region were very large indeed. In this region, for large k, the approximations were either very 
good or very oad. Use of large k approximations should evidently be accompanied by supplementary 
studies at small z. 



9. Additional Comments 

It was mentioned earlier that the plural-series approach can be used to gauge the accuracy of 
representations of the function-fitting type which make use of sums of exponentials, as in eq (A4) of 
appendix A. This requires that the Xj be assigned the values of the specified representation, and that 
the plural-series solution have a{= for n > 0. The latter requirement may or may not be satisfied, 
depending on the norm of the function-fitting representation in the plural-series system. Other 
plural-series representations can have smaller norms, particularly when some of the Aj have 
opposite signs and are large. 

When all the constant factors in expressions such as (A4), in appendix A, have the same sign, 
the simple sum of exponentials tends to be the plural-series solution for all N less than or equal to 
the number of moments given correctly by the approximation. Very likely, with properly stated 
conditions one could prove a theorem to this effect. But when this simple property does not hold, 
one can still modify norm values by adjusting the shape of the error function Z)£ through choice 
of the Aj values. Further, one can modify the moments of the unknown function by addition of 
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moments of known terms, with final subtraction of these functions from the plural-series construc- 
tion. A simple case of this would be transferral of some of the negative terms in a function-fitting 
representation to the other side, with corresponding changes in the moments of the function to be 
constructed. But such a simple prescription cannot be recommended as a general procedure 
because it removes exponential components somewhat indiscriminately from the plural-series 
approximating system. 

Selection of Aj values for the examples of the preceding section was based on the assumption 
that the error function D% should have terms corresponding to an elementary function-fitting rule, 
as discussed in appendix A. Cases with peculiar features, such as large terms which nearly cancel, 
probably require more care in specification of the Aj values. But the quality of the approximations 
usually appears to be surprisingly independent of such assignments. 

10. Appendix A. A Plural-Series Computer Program 

A general-purpose computer program which we have called FFPOL was constructed to perform 
calculations such as those described as examples. The following operations are performed by the 
main constituent components: 

(1) Specification of exponential scale factors Xj, by subroutines designated BCOF and SCALE. 

(2) Solution of the basic Lagrange multiplier problem, using a routine called SEP. 

(3) Computation of the error estimate functions D^ and Z)£, with a routine ESEP, which makes 
use of SEP. 

(4) Evaluation of both approximation and error bound data, with a subroutine designated 
UNESEP. 

More specifically, the following are the computations performed by the main subroutines just 
mentioned: 

The scale factors xj are evaluated as a product 

*j=6/s, (Al) 

where the £ ) are pre-specified according to the rule 

3 = 0MJ=1,2, . . .,/, (A2) 

& = J-\ (A3) 

which ensures that these values cover the unit interval in steps which are basically geometric, 
but with roughly uniform coverage. BCOF evaluated the fj. The overall scale factors is evaluated 
by SCALE so that an expression of the type 

2to/*i +1 )*~ ,/ * J (A4) 

i=i 

can correctly fit (/-hi) moments. This calculation also determines the Aj of (3) and (4), by an 
arbitrary rule such as Aj — \f | . 5 

The procedure followed by SEP in evaluating both the A, values and the aJ n coefficients is a 
straightforward evaluation first of the y\ n matrix elements and then the uu* matrix elements. 



5 The effect of this selection on the quality of the representation is not clear, and not as strong as one might expect. 
According to (41), specification of large Aj values for small xj components would appear to emphasize the small z region in 
the D k functions, for example. The data of tables 1 and 2 were obtained with the function fitting parameters, as in (A4), 
in the combination \fjlx k+1 |. Perhaps a more logical choice for the Aj of (3), (4) might be the combination f z /x k+1 , which 
tends to give D k « \fj \lx k+1 , at least for the dominant terms. 
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The inversion of the latter matrix, followed by the matrix multiplications (16) and (13) completes 
the calculation. 

The error estimate function Z)J is evaluated by ESEP with (25), using (19'). The sums are 
carried to 30-40 terms ordinarily, with evaluation of the t/* (z) mostly by recursion. The recursion 
formula, 

Hn^l)(n^2)[U^ 2 -U^ 1 ]=zW^2(n^l)^n^2k^l)[U^ l -U^ n ] 

- (2n + k) (2n + k - 1) [17* - £/*_J , (A5) 

was derived from the differential equation, using eq (17.16'), both given on p. 744 of reference [2]. 
The recursion approach eventually breaks down, particularly for large values of k; but with suitable 
restrictions, and depending somewhat on the size of the argument z, 30-40 terms gives accurate 
results. ESEP evaluates D k N by subtracting from (D k ) 2 terms of the type P n (z) 2 , as shown in (44). 
Evaluation of the approximation data by UNESEP utilizes (3). The (estimated) error bound 
data is calculated using D§ data from ESEP, with Sat taken to be an adequate approximation to 

Finally, the master routine, FFPOL, first calls BCOF and SCALE. This is followed by an 
optional provision to modify the Xj values just determined to obtain a related set which does not 
have the property of fitting (J + 1) moments, through use of (50). Lastly, UNESEP is called to 
give both approximation and (estimated) error bound data. All calculations have used double 
precision (60 bits). 

1 1 . Appendix B. Convergence of U£ Polynomial Representations 

We summarize and complete here the arguments of reference [4] relating to convergence 
of the t/J series representations. 

(1) First, let us note two properties of these biorthogonal functions: 

U*(z) =O(/i" 1 / 2 ),forz^0, (Bl) 

and, also for z ^ 0, 6 

Note that (Bl) does not usually hold at z = 0, where tf*(0) = 0(/i ( *" 1)/2 ). The U* functions could 
be (but have not been) defined as coefficients of the different powers of u in the expansion (B2). 
This expansion converges uniformly for z # on the unit circle | u | ^ 1, except that, given € > 0, 
we must have | u — 1 1 ^ e. 

(2) Next, let us consider an even function F{z), — °o < z < o°, whose Fourier transform has 
singularities only to the right and far left of the equilateral hyperbola shown in figure Bl, a pair 
of regions which we refer to as "I". The Fourier transform of F(z), which we designate 0(77), 
will be an even function of 77 and will actually have twin singularities symmetrically located about 
the imaginary axis. The Fourier inversion integral is 

1 f io ° 
/r(z)= 2^r7j t dW^W' (B3) 



6 In physical problems, singular behavior at z = will be known and usually treated by a special calculation. The 
remainder at 2 = will often, though not always, vanish. Thus, treatment of this point as a special case does not lead to 
any practical difficulties. Incidentally, eq (40) of reference [4] which is the starting point for the arguments leading to this 
statement on convergence, can be derived by induction without reference to other expressions for U* which include singular 
components at z = 0. This development makes use of (38) of reference [4], and the commutator of the two operators s and 
{d/ds) k in the integrand of (40) of reference [4]. 
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FIGURE Bl. The equilateral hyperbola satisfying |t? 2 — 1 1 — |i7 2 | , which intersects the real axis at ±1/ V2. 

The Roman numerals identify Regions I and II. 

Because of the location of the singularities, we can move the path of integration to a curve C 
completely to the right of the right-hand branch of the hyperbola. 

To apply (B2) to this Fourier integral, we make a change of variables u= (1 — tj~ 2 ) , and 
note that the path of integration C can then be divided into a part d located in the region of uni- 
form convergence and a part Cz corresponding to 1 1 — u \ < e. In the case of C 2 , our change of 
variables gives 

h| > 1/Ve. 

The two curve segments give two integrals; and for C\ it is appropriate to use the series represen- 
ation: 



The exponential factor in the integrand of the last term on the right has a magnitude determined 
by the real part of the exponent, and therefore at least as small as 



exp 



v^j 



For z^O, this factor can be made as small as one pleases; and we choose it so that the second 
term on the right of (B4) can be neglected altogether. 

Because of uniform convergence, the integration and summation can be interchanged in the 
first term on the right of (B4), with the result 



where 



F(z) = 2F«U° n (z)e-U, 



*-hi^W 



(B5) 



(B6) 



We thus find that F(z) has pointwise convergent representations of the type (B5). Also, by writing 
<t>(i7) as a transform over F(z), in (B6), we readily change to the more familiar form 



-„= \" dzij«(z)F(z). 

Jo 
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(B6') 



(3) A parallel argument, also making use of the region of uniform convergence of (B2), but with 
a transformation — w = rj 2 /(l — 7j 2 ) and an inversion integral confined to Region II, leads to the 
series form for the transformation kernels of the £/£ systems: For y, z > 0, 

f^r^o(v7T?)= J {ff*(y)e-W*(*)<r*}, (B7) 

Idydz) ™ £t 

where Ko is the modified Bessel function. For the arguments which follow, the important feature 
of this expression is the fact that for y=z >0, the sum on the right converges to a finite value. 

(4) Let some of the singularities of a function $(17) be located in Region II. Then it may be 
possible to choose a scale factor (3 such that the singularities of O^/jS) all lie in Regions I. If so, 
since 

F(z) = |^ J*™ d V e-^\ie<$>( V lf3) , (B8) 

it is clear that a representation of the type 

F(z)=f i F H U°(zlp)erMli> (B9) 

71=0 

exists. But this will not be the case if any singularities of <$>(r)) lie on the imaginary axis or if any 
point of the imaginary axis is a limit point for singularities of &(r)). 

(5) Preceding arguments applying only to the J7JJ system can be generalized to the U% sys- 
tems, with k > 0. If F(z) is even and is known to have the following trend for small z values: 

F(z)=z*f(z), (BIO) 

then/ (2) may be either even or odd. If the singularities of the Fourier transform <p of/ lie entirely 
in Regions I for some scale factor /3, the argument already given for the k = case leads to point- 
wise convergent series representations of the form 

/(z)=2/*E/*(z/0)e-W, (Bll) 

n = 

with 



^Sj^^V^^ffw^www. 



(BID 



Equation (Bll') makes it clear that/(z) can contain singular terms at 2=0 — Dirac delta functions 
and their derivatives — which do not contribute to the coefficients/^ because they are overriden by 
the factor z k . 

(6) Consider a function F(z) with Fourier transform having singularities confined to Regions I. 
This same property will apply also to/(z), if F(z) has the form (B10). If we write expressions for 
finite approximations to/(z), 



Mz)=2ttU k n(z)e-^ (B12) 

n = 

and apply Schwarz's inequality to the difference function/— /v, we obtain the result 

i/-M 2 ^{i (ft) 2 } f {ffS(*)«H«l}«. (B13) 

l n =0 J n=N 



n = N 
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Since the last factor on the right tends to vanish in the limit /V— » oo 9 it is clear that finite norm, 

f (/i') 2 <^ 



n = 



implies pointwise convergence, forz=N0. 

(7) Let F(z) be even and representable by a f/JJ series, 

F(z)="2FnUl(z)e-\z\. (B15) 

n = 

We further assume that this series has finite norm, and calculate the Fourier transform of both 
sides of (B15), by a limiting process: 



<J>(i7) = lim I dze^F (z) -\- \ dze^F(z) 

€->0 [ J-00 h 

= lim V F n \ \ ' Aevl/;WrW+ f°° <fee"*£/J!(z)e-w] 

€-*0 n=o I J-x A J 

00 2 / T) 2 



7J 2 \ T7 2 — 1 



(B16) 



The limit e~ >0 can be taken as above providing F(z) is sufficiently well-behaved nearz=0 that 
the Fourier transform exists. 

Now, applying Schwarz's inequality to (B16) we obtain 

O f x ~2 2n+2^|l/2 

i*wi«iviii{££l } < B "> 

And since |t7 2 /(t7 2 — 1) | < 1, for 77 in Region II, the right side is finite at all points r) in Region 
II. We conclude that for series representations with finite norm, the Fourier transform can have 
no divergent singularities in Region II. 

(8) Note that for z/0, the cases — °° < \z\ < °° and 0< z < °° can be considered equivalent, 
for functions Fd-zl) and F(z), respectively. Hence problems in which the functions are defined 
on the infinite interval, can be analyzed with the use of [/* representations defined on the semi- 
infinite interval. 

(9) Although we do not attempt it here, the existence of simple generating functions for the 
£/£ polynomials makes a similar set of arguments possible, but in connection with even or odd 
functions whose Laplace transforms have singularities confined to Region II of figure Bl. 

(10) Physical solutions to neutron and gamma ray transport problems are such that the singu- 
larities of the Fourier transform lie on the real axis at distances from the origin |tj| ^ Mmin, where 
)Lt min is the attenuation coefficient of the most penetrating radiation component and has a known 
value in any given problem. Thus one always expects convergent representations of the general 
type (Bll) to exist. 

One should realize that there exist functions of \z\ on the infinite interval whose moments all 
vanish. Such functions cannot be represented by C/J expansions; and one infers that these functions 
have Fourier transforms which are singular on the imaginary axis. This observation is confirmed by 
the discussion of these functions in reference [10]. Since there is neither physical nor mathe- 
matical reason to expect that such functions play any role in transport theory, their automatic 
exclusion from any representation based on U k n functions is appropriate. 
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